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Summary 

A method has been presented for aerodynamic de- 
sign for a specified pressure distribution using anal- 
ysis codes only. The method requires a very con- 
servative number of analysis runs, and, therefore, is 
appropriate when the analysis code is a large code in 
terms of storage and/or running time. Three model 
problems illustrate some capabilities and limitations 
of the method. 

Introduction 

For some of the simpler aerodynamics design 
problems (especially two-dimensional), an inverse 
method exists so that a geometry can be directly 
generated from a specified pressure distribution. For 
somewhat more complex problems, an approximate 
inverse can be combined with an “exact” (in some 
sense) analysis code to yield, by iteration, a design 
to the accuracy of the analysis codes (refs. 1 and 2). 
For more complex problems, even an approximate 
inverse may not be available, or, if available, may 
be difficult to implement. Consequently, considera- 
tion has been given to design by the use of analysis 
codes only (refs. 3, 4, and 5). In this case, one ap- 
proach is to treat the design problem as an optimiza- 
tion problem, with the objective function being the 
mean square error between the output pressure and 
the desired pressure. 

Some problems inherent in the application of a 
numerical optimizer to such design problems have 
been discussed in some detail in references 4 and 5. 
The problem that is addressed here, however, is the 
one that arises when the analysis code is large in 
storage and/or running time. Then the many calcu- 
lations required by a numerical optimizer to compute 
local sensitivities and gradually march toward the op- 
timum become prohibitively lengthy and expensive. 
Furthermore, some problems incorporating inequal- 
ity constraints can yield a design that represents a 
mathematical optimum but not an aerodynamically 
improved design. (See discussion in ref. 5.) 

To avoid these problems, full advantage should be 
taken of such things as the designer’s knowledge and 
experience, available experimental data, and previ- 
ous calculations (refs. 4 and 5). One way that this 
information can be utilized is to generate a set of 
candidate shapes that are individually practical, al- 
though not optimal designs, so that an optimal com- 
bination of these shapes is obtained (ref. 3). This 
approach has at least two advantages. First, if the 
geometric constraints are linear and if each candidate 
shape satisfies the constraints, then a linear combina- 
tion of the shapes satisfies the constraints. Therefore, 
the optimization procedure will not produce an aero- 


dynamically impractical design. The second advan- 
tage to this approach is that, if each of the candidate 
shapes is a “good” shape (its performance is not too 
far from optimal), then presumably the optimizer has 
only a short distance to go to obtain the optimum; 
consequently, relatively few calculations are required. 

These shape functions utilized by the optimizer 
are determined by a different process in the method 
of reference 5. There, they are selected to be “aerody- 
namically specific” so that one shape function is cho- 
sen to control airfoil thickness distribution, another 
to control camber distribution, etc. In the present 
procedure, the shapes utilized are all presumed to be 
intelligent guesses at an optimum design. 

One way of representing the design problem is 
that of obtaining the maximum value for the negative 
of the mean square deviation of the obtained pres- 
sure distribution from the desired pressure distribu- 
tion. With this characterization, the present method 
assumes that shape functions generated by the de- 
signer are rough approximations to the maximum, so 
that a sensitivity-type analysis can be applied locally 
near the maximum. In comparison, conventional op- 
timizers only require that the initial guess functions 
lie somewhere on the “hillside” associated with the 
local maximum. 

If the sensitivity analysis is valid near the opti- 
mum, and if the optimum geometry can be closely 
synthesized as a linear combination of the shape func- 
tions, then this optimum combination is generated by 
the analysis. Failure of either of these conditions is 
indicated by the results. In either case, additional 
shape functions must be added to the data base. 
These various situations are illustrated by the three 
model problems presented as examples. 

The local sensitivity analysis near the maximum 
is analytical and uses the concept of base and cali- 
bration (or comparison) functions described in refer- 
ence 6 and applied in references 7 and 8. 

Symbols 

a, a, b, c coefficients used in various series of 
developments 

C p pressure coefficient 

d difference function, Vd~P ( see e Q s - (5)) 

d r = Pd~Pb (see eqs. (5)) 

€{ mutually orthogonal functions 

e 2 - mutually orthonormal functions 

E error function (eq. (8)) 

(F\, F 2 ) inner product of functions F\ and F 2 , 
defined by equation (Al) 


/ function denoting input to computer 

program 

L linear operator defined by equation (1) 

M free-stream Mach number 

p static pressure distribution 

Poo free-stream pressure 

Q[f] nonlinear operator symbolically 

denoting output of analysis code 

r radial coordinate 

x , y, 2 Cartesian coordinates 

y value of y normalized with respect to 

edge value 

a incidence angle 

Subscripts: 

6 base function 

d desired, or design, pressure distribution 

n denotes quantity in nth term of series 

Analysis 

As in reference 8, the analysis code, or system 
of codes, is treated as an operator operating on the 
input function. Furthermore, the linearity assump- 
tion for input increments is made. This assumption, 
which is the basis for sensitivity analysis, requires 
that the variations in the output of the program be 
approximately linear functions of the corresponding 
small increments in the input function. When this 
condition is not met, as, for example, in problems 
involving shock waves, the problem can sometimes 
be preconditioned so that it becomes linear in small 
variations. Such a preconditioning technique by co- 
ordinate stretching is described in reference 7. 

The linearity condition can be expressed symbol- 
ically (as in ref. 8). Let / denote the input function 
and let Q[f ] denote symbolically the nonlinear oper- 
ator represented by the computer code. Then, it is 
assumed that Q behaves as a (Gateaux) differentiable 
operator, so that, for increments 6f in /, the vari- 
ation in the operator is approximately represented 
by its first variation, L, which is linear in 6f. In 
equation form, 

Q[fo + 6f]-Q[fo] = L[Sf] ( 1 ) 

In other words, for slightly different input functions 
f () and f 0 + 6f, the difference in the output is linear 
in 6f. 


The method to be described is fairly general 
in concept. However, for purposes of simplicity 
in presentation it will be discussed in terms of a 
typical problem — that of designing a geometry to 
obtain a close approximation to a desired pressure 
distribution. 

Assume then that several configurations, repre- 
sented by the inputs ( i = 1 , . . . , n) have been 

analyzed so that the corresponding output pressure 
distributions p 2 - are available. It is assumed that the 
fi functions are aerofunction shapes (in the sense of 
ref. 4), that each satisfies the required geometric 
constraints, and that each f t has been initially se- 
lected so that the corresponding pi roughly approxi- 
mates the desired pressure distribution p^. 

One might attempt to represent the desired dis- 
tribution Pd as a linear combination of the p^, 

n 

Pd^Y “iPi ( 2 ) 

l 

However, some constraint is normally applied to the 
d{ as a result of geometric constraints on the input 
geometries. The most natural constraint, as is shown 
later, is to require 

n 

Y 5 * = 1 ( 3 ) 

l 

This constraint reduces the number of arbitrary co- 
efficients in equation (2) to n— 1. Thus, a new set of 
n — 1 coefficients a ^ (i = 1, . . . , n— 1) is defined by 

n— l 

Pd~Pn~Y a i (P n ~ Pi> ( 4a ) 

1 

This equation can also be written in the form 

n— 1 

Pd^Pn+Y a i ( P n ~ Pi) ( 4b ) 

1 

or 

( n — 1 \ n — 1 

1 + Y a iJ P n ~ Y a iP* ( 4c ) 

Comparing equation (4c) with equation (2) yields 
the relation between the d{ and a z - coefficients. (See 
eq. (15).) In the form of equation (4c), it is clear that 
the coefficients sum to one. 

In equation (4a), p n is denoted the base function 
and the remaining p^ values are the calibration, or 
comparison, functions. Actually, any one of the 
computed p t can be selected as the base function and 


2 


I 


r 


the remaining functions renumbered from 1 to n— 1. 
Thus, for the representation 


n — 1 

Pd~Pb*Y a i (Pt> ~ Pi) (5a) 

1 

or with the difference notation 

n— 1 

d r ^Y (5b) 

1 

the problem is to determine the a z so that the right- 
hand side is the closest possible approximation to 
the left-hand side. If the representation functions 
d t were orthogonal, the a { would be determined by 
the Fourier coefficient formulas in order to yield the 
minimum square error. In general, these functions 
are not orthogonal, but they can be orthogonalized 
and normalized by the well-known Gram-Schmidt 
procedure (ref. 9, p. 116, for example). That is, a 
set of orthonormal functions e 2 - (i = 1, . . . , n— 1) can 
be determined as linear combinations of the 

i 

e i — ^ v c ijdj (6) 

i= i 

The details of calculating the c ZJ by the Gram- 
Schmidt method are straightforward but rather 
lengthy and therefore are relegated to the appendix. 

Now d r is approximated by a linear combination 
of the orthonormal functions 


n— 1 

d r ^ ^ ^ (7) 

1 

Then, the coefficients b z are determined by the 
Fourier formula in order to obtain the best approx- 
imation, in the mean square sense, to d r . That is, 
with the terminology defined by equations (Al) and 
(A2) for the error function 


n— 1 

E = d r -Y (8) 

1 

the quantity \E\ is minimized when the b z coefficients 
are obtained from the formulas 

= (dr , ) (9) 

To be useful, the approximation in equation (7) must 
be expressed in terms of the original differences 


this is accomplished by using equation (6): 


n— 1 


n— 1 i 


*y ^ €i — y ^ y ^ Cikdk 

i i= i /t=i 

n— 1 n— 1 

= Y, b * 53 ° ikdk 

t=i fc=i 

n— 1 /n- 1 \ 

= YL ( biCik ) dk 


with 


fc=i \i=i 
n - 1 

= Y2 akdk 

k = 1 


a k ~ Y'i ^ i c ik 
2—1 


( Cik = 0 for k > i) 
(cifc = 0 for k > i) 
(10a) 

(10b) 


Now, replacing the differences in equation (7) with 
the pressure functions (eqs. (5b) and (5c)) and sub- 
stituting from equations (10) into the result yield 


n— 1 

Pd~Pb^Y a k(Pk-Pk) ( 11 ) 

1 

with afc given by equation (10b). If the pressures in 
equation (11) are written as operators operating on 
the shape functions, it has the form 


n— 1 

Q (fd) - Q ( fb ) *Y a klQ ( fb ) - Q ( f k )] ( 12 a) 

1 

Applying the linearity assumption (eq. (1)) yields 

n— 1 

L(fd~fb)*Y a kL(fb-fk) 


= L 


n— 1 


Yl a k ifb fk) 


L 1 


(12b) 


If it is assumed that, locally, L is uniquely invertible 
then 

71—1 


or 


fd- fb™ Y2 a k{fb~fk) 

l 

( 71—1 \ 71—1 

1 + a kJ fb~ Y a kfk 

Equation (13b) can also be written as 

71 

fd ~ Y h kfk 


(13a) 

(13b) 

(14) 
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where 

r ~ a k (A: = 1 , . . . , n— 1) ' 

fit. = n— 1 (15) 

1+ E {k = n) 

i 

Equations (14) and (15) yield the optimum com- 
bination of the shapes /j, i = 1, . . . , n subject to the 
constraint of equation (3). The significance of this 
constraint is now apparent. It constrains the overall 
size of the combination shape 

n 

Yl^kfk 

1 

If, for example, the shapes were axisymmetric fore- 
body shapes with a fixed base radius, the optimum 
combination would also have the required base ra- 
dius. Similarly, if the slope at the nose were con- 
strained to have a fixed value so that each /*. had this 
required slope, then the optimum combination would 
have the required slope. Such constraints, though ba- 
sic, are nonlinear, and so are not automatically sat- 
isfied by a linear combination of shapes that individ- 
ually satisfy the constraints. Hence, the requirement 
for equation (3). 

Equations (14) and (15) give the optimum combi- 
nation of the original shapes — a shape that will yield 
the closest approximation to the desired pressure dis- 
tribution subject to the constraint of equation (3) 
and within the accuracy of the sensitivity analysis 
assumption of local linearity. A computer program 
that calculates the coefficients a ^ requires, at most, 
a few seconds running time on a modern computer. 
In comparison, the analysis program that computes 
the pressures may require minutes to perhaps hours 
if, for example, a Navier-Stokes code is used. 

After the pressure has been calculated for the op- 
timum combination a better approximation can be 
obtained by adding new geometries to the data base. 
The relative values of the coefficients defining the op- 
timum combination shape should provide some in- 
sight as to the manner in which the geometry should 
be varied in order to improve the approximation. 

The mathematics of this derivation is exact. The 
extent to which the optimal shape derived actually 
approximates the desired result depends primarily 
on two factors. The first of these is the extent 
to which the original data base provides the types 
of shapes that are effective in approximating the 
required result, subject to the applied constraints. 
The influence of this factor is seen in the results 
obtained for the second model problem discussed in 
the following section. 


The second factor depends on the accuracy of 
the local linearity assumption. If the difference 
functions involved in the derivation of the coefficients 
do not have small amplitudes, the linearity condition 
is violated. Local nonlinearity may also be associated 
with the manner in which constraints are (or are not) 
applied. The influence of this factor is seen in the 
results obtained for the first model problem discussed 
in the following section. 

Results on Model Problems 

Three model problems were selected to illus- 
trate the application of the theory. The calcula- 
tions involved are much simpler than the type of 
complex calculation for which the method was in- 
tended. However, simplicity in model problems per- 
mits one to concentrate on fundamental concepts 
without lengthy explanation of extraneous consider- 
ations. Thus, although the power of the technique 
is not adequately demonstrated, the problems se- 
lected are sufficient to illustrate the application of 
the method, its accuracy, and some problems that 
can arise. 

The first problem is that of designing a body for 
a specified upper-surface spanwise pressure distribu- 
tion as a function of cross-section shape. The design 
conditions were M = 1.4, a = 5°, body planform de- 
scribed by the equation y = 0.5 x(l— x), Length = 10 
units, pressure specified at x = 7.5, and body semi- 
height to span ratio at x = 7.5 is 1/3. Each of four 
data base calculations satisfied these conditions, but 
with cross-section shapes specified by the equations 



where y is the y coordinate normalized with respect 
to its edge value. They are shown in figure 1(a). The 
pressures were computed by using the finite-volume 
Euler code described in reference 10. Figure 1(b) 
shows the design pressure distribution that was spec- 
ified and the pressure distributed for the four base 
shapes. Figure 1(c) shows the design pressure dis- 
tribution, the theoretically predicted best approxi- 
mation that can be obtained with the four available 
results (from eq. (2)), as well as the actual result ob- 
tained by computing the pressure for the optimum 
combination shape. This cross-section shape is also 
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shown in figure 1(c). Figure 1(d) compares the de- 
sign distribution with that obtained for the optimum 
combination and with the closest approximation ob- 
tained with a single data base shape. 

The only geometric constraint that was applied 
explicitly to the cross-section shape was the require- 
ment that the body semi-height be 1/3 the span at 
x = 7.5. Since each of the data base shapes satisfies 
this condition, the optimum combination shape also 
satisfies the condition. 

However, no geometric constraint was applied at 
the cross-section edge. The four data base shapes 
represent a considerable range of edge curvatures. 
This factor contributes to the considerable nonlinear- 
ity occurring in the calculations, as is demonstrated 
by the significant difference between the actual and 
predicted pressures for the combination shape. The 
edge represents the leading edge for the crossflow, 
and, consequently, its curvature has a dominating ef- 
fect on the spanwise pressure distribution. In fact, if 
the optimal shape had a relatively small edge radius, 
a crossflow shock would be generated; in which case, 
the linearity would be completely destroyed. A small 
edge radius would also present the problem of upper 
surface separation. The shapes analyzed avoid this 
problem but just marginally, a factor that may also 
contribute to the nonlinearity. However, even with 
the error that results from the nonlinearity, the cal- 
culation provides useful information, and the result- 
ing pressure distribution represents an improved ap- 
proximation to the specified distribution. The shape 
generated could now be used as a base shape with 
new comparison functions generated as small pertur- 
bations of this shape. 

The second model problem was to design a body 
of revolution for a low-speed pressure distribution 
specified at a = 0° over the forward 60 percent 
of the body. Four shapes, shown in figure 2(a) 
were analyzed by a panel method (ref. 11), with 
the pressure distribution for each data base body 
shape shown in figure 2(b). The design pres- 
sure distribution (fig. 2(c)) was specified up to the 
70-percent station only, with the anticipation that 
a smooth fairing beyond the 70-percent station 
would not significantly affect the pressures up to the 
60-percent station. 

The predicted and computed pressure distribu- 
tions for the optimal combination are also shown 
in figure 2(c), along with the optimal shape. It 
is seen that the predicted distribution corresponds 
fairly closely with the computed distribution up to 
the 60-percent station. The difference between the 
specified pressure and the predicted pressure results 
from the failure of the data base shapes to provide 
sufficient variation in their corresponding pressure 


distributions. This difference could be decreased by 
judiciously adding to the data base. 

The third model problem studied was to design 
an axisymmetric supersonic forebody for a specified 
meridian pressure distribution at M = 3.0. Geomet- 
ric constraints were Length-to-diameter ratio = 2.5, 
Semicone angle nose = 16.7°, and zero slope at base 
(x = 1.0). Four configurations were analyzed, with 
shapes described by the equations 

ri (x) = 0.3x — O.lx 3 (0 < x < 1) 

r 2 (x) = ri(x) + 0.04(sm7rx) 2 
7*3 (x) = 7*1 (x) -b 0.1x 3 (l — x) 2 
7\j (x) = 7*1 (x) + 0.1X 2 (1 — X ) 3 

The corresponding pressure distributions, com- 
puted by the method of reference 12 are shown in 
figure 3(a). The specified design pressure distribu- 
tion and the predicted and computed pressure distri- 
butions are all shown in figure 3(b), as well as the 
corresponding body shape. 

Although each of the distributions in figure 3(a) 
falls generally above the specified design distribution, 
the result for the optimum combination matches the 
required distribution closely. 

Furthermore, the fact that the predicted result 
closely approximates the computed result for the 
optimum combination indicates that the problem, 
as formulated, has a high degree of linearity. This 
results from the closely controlled end constraints as 
well as the basic nature of the problem. 

Concluding Remarks 

A method has been presented for aerodynamic 
design using an analysis code only. The method is 
appropriate when the analysis code is a large code in 
terms of storage and/or running time so that mini- 
mizing the number of analysis runs is crucial. The 
designer is required to furnish several intelligently 
selected shape functions whose performances repre- 
sent rough approximations to the optimum. Then 
the method yields the optimum combination of the 
shape functions. This combination is a close approx- 
imation to the actual optimum provided that (1) the 
original shape functions were sufficiently close to the 
optimum and (2) there is sufficient variety in the cor- 
responding pressure distributions so that the desired 
distribution can be closely represented as a linear 
combination of the base distributions. Three model 
problems were presented to illustrate the salient fea- 
tures of the method. 

NASA Langley Research Center 
Hampton, VA 23665-5225 
April 28, 1987 


Appendix 

Gram-Schmidt Orthogonalization and 
Series Expansion of Difference Functions 

The inner product of two functions d, and dj is 
defined by the equation 

(di,dj) - J di(x)dj{x)dx (Al) 

where the integral is evaluated over the relevant 
domain. The corresponding norm is 


This procedure is sufficient to derive the e* functions 
from the d{ functions. However, it is advantageous 
to express the e z - as linear combinations of df 


e* 


“ c ijdj 


(AS) 


From equations (A4), 


(^2,ei) 

C21 MMil 

For the general term, we substitute into the summa- 
tion in equation (A5a) from equation (A6): 


i— 1 


j*. 

II 

? 

a- 

^ . 

(A2) 

e^ di 

Then, the orthonormal set e { is defined recursively 
by the following procedure: 

= <k 

II 

(A3) 

= d % 

= do - {d 2 ,e\)e\ 

(A4a) 


h 

62 N 

(A4b) 

= di 

i- 1 

= dj — ^ ^ (d| , 6j) 6j 
J=1 

(A5a) 

where 

e- = — 

12. 1 

(A5b) 

^t— l,i 


J = 1 

i— 1 i— 1 

■EE 

k— 1 j — 1 
i-1 

■E ; 

k— 1 


E 

3 = 1 


fc=l 

t-1 


fc=l 




(Cjfc = 0 for A: > j) 
(c 3fc = 0 for k > j) 

(A7) 


(A8) 

Finally, using equations (A7) and (A8) with equa- 
tion (A5b) in equation (A6) yields 


(k = i) 


c ik ~~ 


3 = 1 


led (* < 0 


(A9) 


The coefficients in equation (A6) are computed as 
follows. From equation (A3), 

Cll = ^I 


where the coefficients are obtained recursively 
from equation (A8). Thus, with c ^ computed by 
relation (A9), equation (A6) is the result required in 
the main text (eq. (6)). 
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(c) Comparison of optimum combination predicted and computed pressures and design pressure distribution 
with corresponding body shape. 


Figure 2. Concluded. 
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